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We analyse the effect of intrinsic fluctuations on the properties of bistable stochastic sys¬ 
tems with time scale separation operating underl quasi-steady state conditions. We first 
formulate a stochastic generalisation of the quasi-steady state approximation based on the 
semi-classical approximation of the partial differential equation for the generating function 
associated with the Chemical Master Equation. Such approximation proceeds by optimis¬ 
ing an action functional whose associated set of Euler-Lagrange (Hamilton) equations pro¬ 
vide the most likely fluctuation path. We show that, under appropriate conditions granting 
time scale separation, the Hamiltonian can be re-scaled so that the set of Hamilton equa¬ 
tions splits up into slow and fast variables, whereby the quasi-steady state approximation 
can be applied. We analyse two particular examples of systems whose mean-field limit has 
been shown to exhibit bi-stability: an enzyme-catalysed system of two mutually-inhibitory 
proteins and a gene regulatory circuit with self-activation. Our theory establishes that 
the number of molecules of the conserved species are order parameters whose variation 
regulates bistable behaviour in the associated systems beyond the predictions of the mean- 
field theory. This prediction is fully confirmed by direct numerical simulations using the 
stochastic simulation algorithm. This result allows us to propose strategies whereby, by 
varying the number of molecules of the three conserved chemical species, cell properties 
associated to bistable behaviour (phenotype, cell-cycle status, etc.) can be controlled. 
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I. INTRODUCTION 


The networks of interacting genes and proteins that are responsible for regulation, signalling 
and response, and which, ultimately, orchestrate cell function, are under the effect of noise-!--. This 
randomness materialises in the form of fluctuations of the number of molecules of the species 
involved, subsequently leading to fluctuations in their activity. Besides external perturbations, 
biochemical reactions can be intrinsically noisy, especially when the number of molecules is very 
low. 

Far from necessarily being a mere disturbance, fluctuations are an essential component of the 
dynamics of cellular regulatory systems which, in many instances, are exploited to improve cell 
function^. For example, randomness has been shown to enhance the ability of cells to adapt 
and increase their fitness in random or variable environments^—. Random noise also serves the 
purpose of assisting cell populations to sustain phenotypic variation by enabling cells to explore 
the phase spaced ~— U 2 . 

One of the mechanisms that allows noise-induced phenotypic variability relays on multi- 
stability 1314 . The basis of this mechanism was first proposed by Kauffman 1 --, who associated 
phenotypes or differentiated states to the stable attractors of the dynamical systems associated to 
gene and protein interaction networks. In the presence of noise, the corresponding phase space 
generates an epigenetic landscape, where cells exposed to the same environment and signalling 
cues coexist in different cellular phenotypes^. 

Multi-stability is also an essential element in the control of cell response and function via sig¬ 
nalling pathways— In particular, bi-stability as a means to generate reliable switching behaviour is 
widely utilised in numerous pathways such as the apoptosis^, cell survival” differentiation^, and 
cell-cycle progression 21,32 pathways. For example, bi-stability is used to regulate such critical cell 
functions such as the transition from quiescence to proliferation through bistable behaviour asso¬ 
ciated with the Rb-E2F switch within the regulatory machinery of the mammalian cell-cycle^—. 

A common theme which appears when trying to model cell regulatory systems is separation of 
time scales, i.e. the presence of multiple processes evolving on widely diverse time scales. When 
noise is ignored and systems are treated in terms of deterministic mean-field descriptions, such 
separation of time scales and the associated slow-fast dynamics are often exploited for several 
forms of model reduction, of which one of the most common is the so-called quasi-steady state 
approximation (QSS A)22. This approximation is ubiquitously used whenever regulatory processes 
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involve enzyme catalysis, which is a central regulation mechanism in cell function^. In this paper, 
we investigate the effects of intrinsic noise on the bi-stability of two particular systems, namely, an 
enzyme-catalysed system of mutual inhibition and a gene regulatory circuit with self-activation. 
The mean-field limit of both these systems has been shown to exhibit bi-stability^2. The aim 
of this paper is to analyse how noise alters the mean-field behaviour associated to these systems 
when they operate under quasi-steady state conditions. 

We note that this work does not concern the subject of noise-induced bifurcations^. Such 
phenomenon has been studied in many situations, including biological systems. An example which 
is closely related to the systems we analyse here is the so-called enzymatic futile cycles. Samoilov 
et al.— have shown that noise associated to the number of enzymes induce bistability. In the 
absence of this source of noise, i.e. in the mean-field limit, the system does not exhibit bistable 
behaviour. The treatment of this phenomena would require to go to higher orders in the WKB 
expansion, which we do not explore here. 

The issue of separation of time scales in stochastic models of enzyme catalysis has been 
addressed using a number of different approaches. Several such analysis have been carried 
out in which the QSSA is directly applied to the master equation by setting the fast reactions 
in partial equilibrium (i.e. the probability distribution corresponding to the fast variables re¬ 
mains unchanged), and letting the rest of the system to evolve according to a reduced stochastic 
dynamical Other approaches have been proposed such as the QSSA to the exact Fokker-Planck 
equation that can be derived from the Poisson representation of the chemical master equation^. 
Approaches based on enumeration techniques have also been formulated^. Furthermore, Thomas 
et al.— have recently formulated a rigorous method to eliminate fast stochastic variables in monos¬ 
table systems using projector operators within the linear noise approximation^. Methods for 
model reduction based on perturbation analysis have been developed in^22. Additionally, driven 
by the need of more efficient numerical methods, there has been much activity regarding the de¬ 
velopment of numerical methods for stochastic systems with multiple time-scales^^. Several of 
these methods are variations of the stochastic simulation algorithm's^, or the r-leap method^ 
where the existence of fast and slow variables is exploited to enhance their performance with re¬ 
spect to the standard algorithms. Another family of such numerical methods is that of the so-called 
hybrid methods, where classical deterministic rate equations or stochastic Langevin equations for 
the fast variables are combined with the classical stochastic simulation algorithm for the slow 
variables^!!. Other related methods were studied in^“—. 
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Here, we advance the formalism developed in^, in which a method based on the semi-classical 
approximation of the Chemical Master Equation allows to evaluate the effects of intrinsic random 
noise under quasi-steady conditions. In our analysis of the Michaelis-Menten model of enzyme 
catalysis in-^, we showed that the semi-classical quasi-steady state approximation reveals that 
the velocity of the enzymatic reaction is modified with respect to the mean-field estimate by a 
quantity which is proportional to the total number of molecules of the (conserved) enzyme. In 
this paper, we extend this formalism to show that, associated to each conserved molecular species, 
the associated (constant) number of molecules is a bifurcation parameter which can drive the 
system into bi-stability beyond the predictions of the mean-field theory. We then proceed to test 
our theoretical results by means of direct numerical simulation of the Chemical Master Equation 
using the stochastic simulation algorithm 54 . We should note the Hamiltonian formalism derived 
from the semi-classical approximation is formulated on a continuum of particles, which requires 
the number of particles to be large enough. This must hold true for all the species in our model, 
both fast and slow. Since this separation between fast and slow species is based on their relative 
abundance, one must be careful that the scaling assumptions are consistent, particularly in the 
case of the model of self-activating gene regulatory circuit where the number of binding sites is 
typically small. This assumption, however, has been used in previous studies^. Also we show that 
our simulation results of the full stochastic processes agree with our analysis and, therefore, our 
re-scaled equations are able to predict the behaviour of the system. We note that the mean-field 
limit, which is conventionally obtained by ignoring noise in the limit of large particle numbers, is 
obtained by setting the momenta in our phase-space formalism to 1. 

The approximation we develop in this paper falls within the general framework of the optimal 
fluctuation path theory^ This framework is a particular case of the large deviation theory which 
allows us to study rare events (i.e. events whose frequency is exponentially small with system 
size). Within these framework we will show that, upon carrying out the QSSA, the only source 
of noise in the system is associated to the random initial conditions of the species whose numbers 
are conserved. We therefore predict that a population of cells, each having a random number of 
conserved molecules, will have a bimodal distribution. 

This paper is organised as follows. Section 2 is devoted to a detailed exposition of the semi- 
classical quasi-steady state approximation for stochastic systems. In Sections 3 and 4, we apply 
this formalism to analyse the behaviour of a bistable enzyme-catalysed system and a gene regula¬ 
tory circuit of auto-activation, respectively. We will show that our semi-classical quasi-steady state 
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theory allows us to study the effect of intrinsic noise on the behaviour of these systems beyond the 
predictions of their mean-field descriptions. We also verify our theoretical predictions by means 
of direct stochastic simulations. Finally in Section 5, we summarise our results and discuss their 
relevance. 

II. SEMI-CLASSICAL QUASI-STEADY STATE APPROXIMATION 

Our aim in this paper is to formulate a stochastic generalisation of the quasi-steady state ap¬ 
proximation for enzyme-catalysed reactions and simple circuits of gene regulation and use such 
approximation to determine if the presence of noise has effects on the behaviour of the system be¬ 
yond the predictions of the corresponding mean-field models. Specifically, we analyse stochastic 
systems for which the mean-field models predicts bi-stability and investigate how such behaviour 
is affected by stochastic effects. Our analysis is carried out in the context of Markovian models 
of the corresponding reaction mechanisms formulated in terms of the so-called chemical master 
equation (CME)^. Two example of such stochastic systems, a bistable enzyme-catalysed system 
and a gene regulatory circuit of auto-activation, are formulated and analysed in detail in Sections 
M and |IVl respectively. Following^, we formulate the QSS approximation for the asymptotic 
solution of the CME obtained by means of large deviations/WKB approximations^—. The CME 
is given: 



( 1 ) 


where W t (X ) is the transition rate corresponding to reaction channel i and r, is a vector whose 
entries denote the change in the number of molecules of each molecular species when reaction 
channel i fires up, i.e. P{X(t + At) = X(t) + r^xit)) = Wi(X)At. 

An alternative way to analyse the dynamics of continuous-time Markov processes on a dis¬ 
crete space of states is to derive an equation for the generating function, G(pi,... ,p n ,t) of the 
corresponding probabilistic density: 



( 2 ) 


X 


where P(X i,..., X n , t) is the solution of the Master Equation ([[])• G(pi,... ,p n ,t) satisfies a 
partial differential equation (PDE) which can be derived from the Master Equation. This PDE is 
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the basic element of the so-called momentum representation of the Master Eq uati o " — . 

Although closed, analytic solutions are rarely available, the PDE for the generating function 
admits a perturbative solution, which is commonly obtained by means of the WKB method^. 
More specifically, the (linear) PDE that governs the evolution of the generating function can be 
written as: 


— = H k (p l ,...,p n ,d Pl ,...,d Pn )G{p 1 ,...,p n ,t) (3) 

where the operator Hk is determined by the reaction rates of the Master Equation ([I])- Furthermore, 
the solution to this equation must satisfy the normalisation condition G(p\ — 1,.. ., p n — 1, t) — 1 
for all t. This PDE, or, equivalently, the operator H, are obtained by multiplying both sides of the 
Master Equation © by H" =1 P?’ an d summing up over all the possible values of (X b ..., X n ) 
From the mathematical point of view, Eq. © is a Schrodinger-like equation and, therefore, 
there is a plethora of methods at our disposal in order to analyse it. In particular, when the fluctu¬ 
ations are (assumed to be) small, it is common to resort to WKB methods^^M. This approach is 
based on the WKB-like Ansatz that G(p\,... ,p n , t) = e~ s< ' pi ’-’ Pn ’ t \ By substituting this Ansatz 
in Eq. © we obtain the following Hamilton-Jacobi equation for the function S(pi ,..., p n , t ): 


dS TT ( dS dS\ 

n, I Pli ■ ■ ■ j Pm Tj j ■ ■ ■ j Tj I 

Ot \ dpi OPnJ 


( 4 ) 


Instead of directly tackling the explicit solution of Eq. ©, we will use the so-called semi- 
classical approximation. We use the Feynman path-integral representation which yields a solution 
to Eq. © of the typ o ^ —— —: 


G(pi,... ,p n , t) = [ e- s ^’- Pn ’ Ql -’ Q ^VQ(s)Vp{s), ( 5 ) 

Jo 

where VQ(s)Vp(s) indicates integration over the space of all possible trajectories and S(p\,..., p n . (),,) 

is given by^: 


S(pi i • • • j Pni Ql i • • • i Qn) J {j>\, j Pm Ql > ■ ■ ■ j Qn) T ^ ^ Qi(s)Pi(s) J ds 

n 

+ ^2s 0 ,i(pi,Qi), ( 6 ) 


2=1 
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where the position operators in the momentum representation have been defined as Q t = d Vi with 
the commutation relation [Qi,Pj\ = So^Sij. S 0ti (pi, Qt) corresponds to the action associated with 
the generating function of the probability distribution function of the initial value of each variable, 
Xi(t = 0), which are assumed to be independent random variables. 

The so-called semi-classical approximation consists of approximating the path integral in Eq. 
©by 


G(pi,. ..,p n ,t) — e~ s ( pi, - ,Pri ti ( 7 ) 

where pi(t),... ,p n (t ) are now the solutions of the Hamilton equations, i.e. the orbits which 
maximise the action S: 

dpi = dH k 
dt dQi 
dQi dH k 

dt dpi 

where the pair ( QuPi ) are the generalised coordinates corresponding to chemical species i = 
1 , ...,n. These equations are (formally) solved with boundary conditions^ Qi( 0 ) = 27(0), 
Pi(t) = Pi, where Xj( 0 ) is the initial number of molecules of species i. 

Eqs. ©-© are the starting point for the formulation of the semi-classical quasi-steady state 
approximation (SCQSSA)^. In order to proceed further, we assume, as per the Briggs-Haldane 
treatment of the Michealis-Menten model for enzyme kinetic s — 1 —, that the species involved in the 
system under scrutiny are divided into two groups according to their characteristic scales. More 
specifically, we have a subset of chemical species whose numbers, X,, scale as: 

Xi = Sxi, ( 10 ) 

where 27 = 0(1), whilst the remaining species are such that their numbers, Xj, scale as: 

Xj = Exj, (11) 

where x 3 =0(1). Key to our approach is the fact that S and E must be such that: 

E 

e = I- (12) 
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We further assume that the generalised coordinates, Qi , scale in the same fashion as the corre¬ 
sponding variable X u i.e. 


Qi = Sq u (13) 

where qi = 0(1). We refer to the variables belonging to this subset as slow variables. Similarly, 

Qj = Eqj, (14) 

where qj = 0(1), which are referred to as fast variables. Moreover, we assume that the moment 
coordinates, Pi, are all independent of S and E, and therefore remain invariant under rescaling. 

Under this scaling for the generalised coordinates, we define the following scale transformation 
for the Hamiltonian in Eq. ©: 

H k (pi, ■ ■ .,p n ,Q l, ■■■,Qn) = k J S k E l H K (p 1 , ...,p n ,q x ,..., q n ) (15) 

where J identifies the reaction with the largest order among all the reactions that compose the 
dynamics and kj is the corresponding rate constant. For example, in the case of the bistable 
enzyme-catalysed system whose reactions or elementary events and the corresponding transition 
rates are given in Table [Q J = 1, as this reaction is order 3 whereas all the others are order 0, 1, or 
2. In the case of the self-activating gene regulatory circuit, Table ITVl J = 3, since this reaction is 
order 3 whereas the remaining ones are order 1 at most. The exponents k and l correspond to the 
number of slow and fast variables involved in the transition rate Wj, respectively. 

The last step is to rescale the time variable so that a dimensionless variable, r, is defined such 
that: 


T = k J S k ~ 1 E l t (16) 

It is now a trivial exercise to check that, upon rescaling, Eqs. ©-© read 

dpi _ dH K 
dr dqi ’ 
dqj _ dH K 
dr dpi 

for the slow variables. By contrast, rescaling of the Hamilton equations corresponding 
subset of fast variables leads to: 


(17) 

(18) 
to the 
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x 2 +x 4 =’— x 6 —x,+x 4 
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0-X 7 - 0 

FIG. 1. Reactions for the bistable enzyme-catalysed system proposed by Tyson & Novak? 2 , X\ repre¬ 
sents active Cdh/Apc, Xi inactive Cdh/Apc, X-> inactivating enzymes, X 4 activating enzymes, AT, active 
Cdh/Apc-inactivating-enzyme complexes, Xq inactive Cdh/Apc-activating-enzyme complexes, and X-j the 
number of CycB-CDK complexes. The first two reactions correspond to enzyme-catalysed inactivation and 
activation of Cdh/APC. The third reaction corresponds to the dynamics of CycB activity: synthesis at a 
constant rate, k-j, and degradation by natural decay and active Cdh/Ape-induced inactivation. 


dpj _ dll K 
dr dqj ’ 
dqj _ dll i; 

6 dr dpj ’ 


(19) 

( 20 ) 


where e is defined in Eq. (fT2l) . The QSS approximation consists on assuming that ~ 0 and 
e ^± ~ 0 in Eqs. ([I9t-(l20b. 


dH K 
dqj 
0H k _ 

dpj 


= 0, 


0, 


( 21 ) 

( 22 ) 


resulting in a differential-algebraic system of equations which provides us with the semi-classical 
quasi-steady state approximation (SCQSSA). 


III. BISTABLE ENZYME-CATALYSED SYSTEMS 

As a prototype of a bistable enzyme-catalysed system, we analyse a stochastic system proposed 
in^iZS, whose mean-field limit has been shown to correspond to a bistable system which is a part 
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Variable Description 


X\, X 2 Number of active and inactive (respectively) Cdhl molecules 
X 3 , X 4 Number of Cdhl-inactivating and Cdhl-activating (respectively) enzyme molecules 
X 3 , X(, Number of enzyme-active Cdhl and enzyme-inactive Cdhl (respectively) complexes 
Xj, Number of active cyclin molecules 


Transition rate 


Event 


W\{x) = k\X' 7 X\X i n = 

W 2 (x) = k 2 X 7 X 5 r 2 = 

W 3 {x) = k 3 X 7 X 5 r 3 = 

W 4 (x) = k 4 X 2 X 4 r 4 = 

W 5 (x) = k 5 X 6 r 5 = 

W 6 (x) = k 6 X 6 r 6 = 

W 7 (x) = k 7 r 7 = 

W 8 (x) = k 8 (l + aX 1 )X 7 r 8 = 


(—1,0, —1, 0, +1, 0,0) Enzyme and active Cdhl form complex 
(+1,0, +1, 0, —1, 0,0) Enzyme-active Cdhl complex splits 
(0, +1, +1, 0, —1, 0,0) Inactivation of Cdhl and enzyme release 
(0, —1,0, —1,0, +1,0) Enzyme and inactive Cdhl form complex 
(0, +1,0, +1,0, —1,0) Enzyme-inactive Cdhl complex splits 
(+1,0,0, +1,0, —1,0) Activation of Cdhl and enzyme release 
(0,0, 0, 0,0,0,+1) CycB synthesis 
(0,0, 0,0,0,0, —1) CycB degradation 


TABLE I. Random variables and transition rates of the stochastic model associated to the enzymatic reaction 
shown in Fig. [U 


of a model for the Gi/S transition of the eukaryote cell cycle proposed in- 22 . Tyson & Novak- 2 have 
formulated a (deterministic) model of the cell cycle such that the core of the system regulating the 
G]/S transition is a system of two mutually-repressing proteins (Cdhl and CycB). This system of 
mutual repression gives rise to a bistable system where one of the stable steady states is identified 
with the Gi phase whereas the other corresponds to a state where the cell is ready to go through the 
other three phases of the cell-cycle, known as S, G 2 , and M. This central module, which is the one 
we focus on, is acted upon by a complex regulatory network which monitors if conditions are met 
for the cell to undergo this transition and accounts for its accurate timing. Presently, we ignore this 
network and focus on the central bistable system. It is shown in=^ that the mean field version of the 
model exhibits bistable behaviour as a function of a bifurcation parameter m, i.e. the mass of the 
cell. For very small values of m, the system is locked into a high (low) Cdhl (CycB)-level stable 
fixed point (i.e. into the Gi phase). For very large values m, the system has only one stable steady 
state corresponding to a low (high) Cdhl (CycB)-level fixed point. For intermediate values of m 
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the system exhibits bistability, i.e. both of these stable fixed points coexist with an unstable saddle 
point. In this section, we focus on how noise alters the behaviour of the mean-field dynamics. 

The transition rates corresponding to the different reactions involved in the stochastic model 
associated to the enzyme-regulated ki netics shown in Fig. |Tj are given in Table HI This kinetics 
corresponds to the enzyme regulated activation and inhibition of Cdhl (an inhibitor of cell-cycle 
progression). Cdhl inactivation is further (up)regulated by the presence of CycB, an activator of 
cell-cycle progression. CycB is synthesised and degraded at basal rates and is further degraded in 
the presence of active Cdhl (see Fig. Q]). Therefore, the resulting dynamics leads to a system with 
mutual inhibition which produces bistable behaviour. It is important to note that the associated 
reaction kinetics exhibits three conservation laws (see Table [D): X 3 + X 5 = e 0 , X 4 + X 6 = e 0 , 
and X t + X 2 + X 5 + X (i = s 0 . The first two of these conservation laws are associated to the 
conservation of the number of Cdhl-inhibiting and Cdhl-activating enzymes, respectively, whilst 
the latter expresses the conservation of the total number of Cdhl molecules. The quantities e 0 
and s 0 are the (conserved) number of enzymes and Cdhl, respectively. Note that, as per the 
methodology developed in Section UH we assume that s 0 = 0(S ) and e 0 = 0(E). 


Rescaled variables 

Dimensionless parameters 

r = k\ESt 

e = E/S, a = aS 

<7i = Qi/S 

k 2 = k 2 /(hS) 

<72 = Q 2 /S 

K 3 = h/(h S) 

Q3 = Q'i/E 

n 4 = k 4 / ( k\S ) 

<74 = Qa/E 

= k 5 /(kiS 2 ) 

Q5 = Qb/E 

^6 = k 6 /(kiS 2 ) 

Qg = Qg/E 

k 7 = kj/ihES 2 ) 

<77 = Qi/S 

K8 = ks/(hES) 


TABLE II. Dimensionless variables used in Eqs. (l29l) . S and E are the average concentration of Cdhl 
(active plus inactive) and the average concentration of both Cdhl-activating and Cdhl-inactivating enzymes, 
respectively. We further assume that the stationary concentration of active CycB also scales with S. 

The corresponding stochastic Hamiltonian, H k , which is derived by applying the methodology 
of Section [TQ to the Master Equation associated to the chemical kinetics described in Table HI can 
be split into three parts, 


12 




Hk(pi, ■ ■ ■ ,P7, Q i, ■ ■ ■, Qi) — Ha + Hi + H B , (23) 

where Hi is the Hamiltonian corresponding to the CycB-regulated enzymatic inactivation of Cdhl 
(reactions 1 to 3 in Table U): 

Hi(p, Q) = h(p 6 - p 2 p^)Q 2 Qi + k 5 (p 2 p 4 ~ Ps)Qq + h{PiPi -Pa)Q6, (24) 

Ha corresponds to enzymatic activation of Cdhl (reactions 4 to 6 in Table U): 


H a (p, Q ) = hp 7 (p 5 - Pip 3 )QiQ 3 Q 7 + k 2 p 7 {pip 3 - P^QsQi + k 3 p 7 (p 2 p 3 - p 3 )Q 3 Q 7 i (25) 

and, finally, H B , which corresponds to synthesis and degradation of CycB, is given by (reactions 
7 and 8 in Table U): 

H B (p,Q ) = k 7 (p 7 - 1) + A*(l ~p 7 )Q 7 + k 8 api(l -p 7 )QiQ 7 . (26) 

We now proceed to apply the procedure explained in Section 2 in order to obtain the SCQSSA 
for the system determined by the transition rates given in Table [0 We first need to determine 
which of the variables are slow variables and which ones are fast variables. As shown in Table 
HH the pairs (pi,Qi), {p 2 ,Q 2 ), and (p?, () 7 ), corresponding to the active and inactive forms of 
Cdhl and to CycB, respectively, are the slow generalised coordinates, as the generalised positions 
scale with s 0 . The remaining generalised coordinates scale as e 0 and, therefore, are fast variables. 
Furthermore, the rescaled Hamiltonian is given by: 

H k (p,Q) = k 1 ES 2 H K (p,q ) (27) 

where 

H K (p, q) = H k A + H^i + H KiB , (28) 

with 


H K J = K 4 (P6 - _P2P4)<?2<?4 + K 5 (p 2 P4 - Pe)?6 + «e (PlPl - Pd)q 6 

H k ,A = Pt(P 5 - PlP 3 )gi<? 3<?7 + K 2 p 7 (pip 3 - £>5)7.577 + K 3 p 7 (p 2 p 3 - £>o) 7 o 77 

H k ,b = K 7 (p 7 - 1) + k 8 (1 - ^7)77 + K 8 api(l - £>7)7177 ( 29 ) 
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The rescaled parameters Ki are given in Table Hfl Last, by rescaling time and defining the dimen¬ 
sionless time variable as r = k\ESt (Table HD, the SCQSSA equations (fT7T)-(fT8l) and (l2TI)-(l22l) 
lead to (see^ for a detailed derivation): 


dqi 

dr 

<kj2 

dr 

<k]_ 

dr 


PiPe A 


« 6 g 2 


«3?79i . n v 

T P7PzPe 3 -—- + K 8 a(l- p 7 )q 7 qi 

<?2 + d 2 q\ + Ji 

n e q 2 , - " " 


= «7 - « 8 (1 + apiqi)q 7 


P 5 = P 3 Pl 


P6 = P4P2 

-r- = “(I - Pr)«s(l + apiqi) 


(30) 

(31) 

(32) 

(33) 

(34) 

(35) 


where p 1 = p 2 , p 3 , and p 4 are constants to be determined and J x = n 2 + k 3 and J 2 = kJ 1 (k 5 + n 6 ), 
andpe 3 = e 3 /E andp e4 = e 4 /E. Note that for qi(r) +q 2 (r) = p c , withp c = s 0 /S, to holdp 7 = 1 
must be satisfied. In this case, we have 


dqi K 6 (p c ~q\) k 3 mq 7 q 1 

dr \Pc~ qi) + J 2 qi + Ji 

(36) 

~r = k 7 - k 8 (1 + oip\qi)q 7 

(IT 

(37) 

P5 = P 3 Pl 

(38) 

Pe = PaPi 

(39) 


As shown in^, the parameter values are determined by comparing the corresponding mean- 
field approximation, which is obtained by taking pi = 1 —, and p c = p e3 = p f i = 1 , i.e. the total 
number of molecules of Cdhl and its activating and inhibiting enzymes be exactly equal to its 
average, i.e. s 0 = S and e 3 = e 4 = E, to the system originally proposed by Tyson & Novak—. In 
Eq. (l36l) we have redefined k , 3 —* n 3 m in order to make explicit the dependence on the bifurcation 
parameter, m, as used by Tyson & Novak—. The parameter values are shown in Table HID 

Upon rescaling of the variables (Table HD and the Hamiltonian (Eq. ([131) 1. the action functional 
reads: 
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(a) 




(b) 



FIG. 2. (a) Bifurcation analysis for the SCQSS approximation of the stochastic bistable enzyme-catalysed 
system Eqs. (|36T)-(|391). The panels on the top plot (a) shows the bifurcation diagrams for different values 
of the parameters pi, p c = p\ and p = 3 . If eo and so are random Poisson variables with parameter S 

and E, respectively, then p = -§ (see Eq. 1451) . In these panels solid lines correspond to r = 1, dot-dashed 

P-i 

lines to r = 2, and dashed lines to r = 3. The bottom plot (b) shows the bi-stability boundaries in p\ — m n 
parameter space. The region between the boundaries corresponds to the bistable region of the stochastic 
Tyson & Novak system according to the SCQSS approximation. 
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Rescaled parameter Parameter 


Units 


Reference 


k 2 = J4 — k 3 

= a ± m 
k\ES 


Ct/Q 

K 6 = - 1 


a 3 


K5 — H4J3 — ^6 


u -1 

= — 


fciSS 


U/o 

K S = — 


fclES 


u , 2 

° = fci £Sk 8 


a; = 0.04 


a' 2 = 0.04 
a 2 = 1 


a 3 = 1 
a4 = 35 
m = 0.3 
E = 0.01 
5 = 1.0 
fci = 1 
k 4 = k 3 
J 3 = J 4 = 0.04 


mi n 


mm 


mi n 


mm 


-1 

-1 

-1 

mm -1 
Dimensionless 
Dimensionless 
Dimensionless 
mm -1 
Dimensionless 
Dimensionless 


TABLE III. Parameter values used in simulations of the stochastic bistable enzyme-catalysed system 


S(p, q )=s 0 j 

U \ slow fast 

n 

+ y ] Sp ,i(pi) (40) 

i 

It is straightforward to check that in SCQSSA conditions H K ( 79 . q) = 0. Furthermore, since p\ = 
p 2 =const. and p 7 = l,and epj ~ 0 for the fast generalised coordinates, the SCQSS approximation 
of the action Eq. (1401) . Sq S s , reduces to: 

n 

Sqss(p) = ^2s 0 ,i(pi) (41) 

i =1 

where, as per the SCQSSA, p 5 and p 6 are determined by Eqs. (l38l) and ([39b . respectively, p- = 1, 

which implies S 0}7 (p 7 ) = 0, and p\ = p 2 , P:>, and p 4 are constants that remain to be determined. In 

order to do so, we resort to the method developed in reference^. The quasi-steady state character¬ 
istic function, Gqss(p , t) is given by: 

6 

G Q ss(p,r) = = J]G 0 .(p.) (42) 

1=1 
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where G 0)i (pi) = e~ s °’ i ( pi ' ) is the generating function of the probability distribution for the initial 
condition of species X t i = 1,..., 6. In^, we have shown that, applying a Laplace-type asymptotic 
methodZLZS to the integrals 



£ 


e -(So,i(Pi)+ s o log pi) 


dp i, 


Pi 



1 f e _ ^ So ’^ p ^ +e ° logp ^ 


£ 


dpi with i — 3,4, 


(43) 


Pi 


where, p 4 = p 2 , P 3 and p 4 can be given as functions of s 0 and e it i — 3,4, i.e. the initial numbers 
of Cdhl molecules and Cdhl-inactivating and Cdhl-activating enzymes, respectively: 


dS Qd 



(44) 


P(Xi(t = 0) = s 0 ), P(X 3 (t = 0) = e 3 ) and P(X 4 (r = 0) = e 4 ) are the probabilities that X 4 
initially takes the value A", (r = 0) = s 0 and that X 3 and X 4 have initial values X 3 (t = 0) = e 3 
and X 4 (t = 0) = e 4 . These probabilities can be interpreted to correspond to variability in the 
abundance of these enzymes within a population of cells. A particularly simple case results from 
assuming that P(X 4 (r = 0) = s 0 ), P(X 3 (r = 0) = e 3 ) and P(X 4 (r = 0) = e 0 ) are Poisson 
distributions with parameter S and E, respectively. In this case^: 


Note that, in the particular case in which the total numbers of Cdhl and enzyme molecules are 
random Poisson variables, we have that p\ = p c . p : > = p e3 , and p 4 = p e4 . 

A. Bifurcation analysis 

Fig. [2] shows results regarding the bifurcation behaviour of the SCQSS approximation of the 
stochastic bistable enzyme-catalysed system Eqs. (I36l)-(l39l). In particular we are interested in 
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FIG. 3. Simulation results for the stochastic bistable enzyme-catalysed system Table U We have plotted the 
probability P{x\,T ) = Prob(xi(r = T)) where x\ = X\/S and T = 100 for different values of p. The 
initial number of Cdhl-inactivating and Cdhl-activating enzymes are fixed according to X 3 (t = 0) = ^ 
and X 4 (t = 0) = eo, respectively, rn = 0.3. We aim to check our predictions regarding the effect of the 
ratio p = -f = -f on the stability properties of the system. According to our results shown in Fig. [2J 

Pi e 4 

decreasing the ratio between the number of Cdhl-inactivating (e 4 ) and Cdhl-activating ( 63 ) enzymes, the 
system should be driven away from bistability and into the stable Gi-phase regime (see Fig. [2jb)). The 
remaining parameter values are inferred from those given by Tyson & Novak 22 as shown in Tables HO and 
m We see that when varying p, the system switches from a state of high x\ (p > 0.9) ro a state of low 
x] (p < 0.6), whereas at the intermediate levels of (e.g. p = 0.7 and p = 0.8) the system is in a bistable 
state. We take pi = p c = 1 in all the simulations shown in this figure. Average is performed over 1000 
realisations. 


a comparison between the bistable behaviour of the mean-field model, corresponding to taking 
Pi = 1 for all i, and that of the SCQSS approximation with p 3 , p 3 and p 4 given by Eq. (l44l) . i.e. 
they are determined as functions of so and eo- 

2 2 

We have shown that both the ratio of p 3 and » 4 , p = ’XX-X = % = %, and p\ alter the 
bistable behaviour of the system beyond the predictions of the mean-field model. In particular, 
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we observe that decreasing the value of p extends the region of stability of the Gi-fixed point, 
i.e. the fixed point corresponding to the steady-state value of q\, such that q 3 ~ 1. By contrast, 
when p is increased the stability region of the Gi-fixed point shrinks. Intuitively, given the relation 
between p 3 and p 4 and the number of Cdhl-inactivating and Cdhl-activating enzyme, this result is 
straightforward to interpret: decreasing the number of Cdhl-inactivating enzyme demands a larger 
value of m in order to de-stabilise the Gi-fixed point. This is fully confirmed by direct simulation 
using Gillespie stochastic simulation algorithm 54 . Fig. [3] shows simulation results in which we 
compute the probability P{x i, T) = Prob(xi(r = T)) for different values of p < 1. T has been 
chosen so that the system has reached steady state conditions. We observe, that for p = 1 and 
m = 0.3, the system evolves towards the q 1 <C 1-fixed point (i.e. the S-G 2 -M fixed point). As p 
decreases, i.e. there is more Cdhl-inactivating enzyme than Cdhl-activating enzyme, the system 
enters the bistable regime. If p reaches low-enough values (depending upon the initial condition), 
we may even observe an exchange of stability, i.e. the system evolves towards the q x ~ 1 -fixed 
point. 

Regarding the dependence on p } , we have checked the predictions of the SCQSS approximation 
by means of simulations with different values of so- Figure [2] shows the bi-stability region of 
system Eqs. (I36l)-(l39l) in p 3 — m R - space, where m R = pm. For a fixed value of m R , there 
is a threshold value for p 3 below which the system stops being bistable to become entrapped 
into the the S-G 2 -M fixed point (i.e. <p <C 1). In order to validate this prediction, we have 
conducted stochastic simulations for different values of so- Figure 0 shows simulation results for 
P(xi,T) = Prob(xi(r = T)). We observe that for small values of s 0 , the system is locked into the 
the S-G 2 -M fixed point, as predicted by the SCQSS approximation. As s 0 increases, the system 
enters a fluctuation-dominated bistable regime where, as the system goes through the bifurcation 
point, the system undergoes bistable behaviour. This behaviour is typical in a system undergoing 
a phase transition, where fluctuations unboundedly increase^. Finally, as s 0 continues to increase, 
the system becomes trapped into Gi-fixed point (see Figure 0). These results fully reproduce the 
behaviour predicted by our SCQSSA stability analysis. 

The aforementioned behaviour regarding unbounded increase of fluctuations close to a bifurcation^ 
is used to locate the critical value of the associated control parameter, i.e. p and p c for the sim¬ 
ulations shown in Figs. 0 and 0 respectively. This property allows us to do a quantitative com¬ 
parison between the simulations and asymptotic analysis. To this end, we plot how the variance, 
a 2 = ((Xi — (xi)) 2 ) where x 3 = X\{r = T)/S, changes as the corresponding control parameter 
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FIG. 4. Simulation results for the stochastic bistable enzyme-catalysed system Table U We have plotted the 
probability P(x\,T) = Prob(x'i (r = T)) where x\ = X\/ S and T = 100 with different initial conditions 
and different values of p c . Average is performed over 1000 realisations, m = 0.3 and X : >(t = 0) = eo 
and X 4 (t = 0) = eo- The remaining parameter values are inferred from those given by Tyson & Novak- 22 
as shown in Tables [FT] and [HI] We see that when varying p c , the system switches from a state of high x\ 
( p c > 0.8) ro a state of low x \ (p c < 0.6), whereas at the intermediate levels of p c = 0.7 the system is in a 
bistable state. 


varies. Regarding the results shown in Fig. Oa) (associated to the simulations shown in Fig. [3]), 
we observe that the critical value of the control parameter p, p B , is approximately p B — 0.7, 
which, taking into account that m = 0.3, implies that the critical value of the renormalized mass, 
m R = pm, m B = p B m ~ 0.21. Our asymptotic analysis predicts that m B = 0.11 (see Fig. [2^b) 
with p c = 1). The results shown in in Fig. [5Jb) (corresponding to the simulations shown in Fig. 0]), 
the critical value of p c , p B , is approximately p B ~ 0.7. The prediction of our asymptotic analysis 
(see Fig. [2];b) with p = 1) is p B = 0.6. 
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(a) 


(b) 




FIG. 5. Plots showing the variance a 2 = ((x'i — (xi)) 2 ) where x\ = X \(r = T)/S associated to the 
simulation results shown in Fig. [3] (panel (a)) and in Fig. 0] (panel (b)). These plots show how <r 2 changes 
as the control parameter ( p , for the simulations associated to plot (a), and p c for the simulations shown in 
plot (b)). The maximum of a 2 as a function of the control parameter helps us to quantitatively determine 
the corresponding critical value 73 . 

IV. AUTO-ACTIVATION GENE REGULATORY CIRCUIT 

We now proceed to analyse the effects of intrinsic noise in a model of a bistable self-activation 
gene regulatory circui t 3 -'—' 74 in the context of the quasi-steady regime. Many instances of genetic 
switches, i.e. bistable gene regulatory circuits, have been identified^ — — 77 . Most of them are 
characterised by the presence of a positive feed-back in which one of the molecular species in¬ 
volved in the system up-regulates its own production. All of these systems exhibit bi-stability and 
hysteresis, i.e. a form of memory associated to bistable systems, and some of them are thought 
to exist in regimes where stochastic switching is frcqucnt^-^. Noise effects on this kind of sys¬ 
tem has been extensively analysed and found to have both constructive and deleterious effects. 
For example, Frigola et al.— have found that noise stabilises the inactive (OFF) steady-state of a 
model of a bistable self-activation gene regulatory circuit by extending its stability region. In this 
Section, we analyse the effects of noise specifically associated to the quasi-steady state regime in 
the large-deviations (large number of molecules) limit. 

We study the stochastic system of the simple self-activating gene regulatory circuit schemati- 
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FIG. 6. Schematic representation of the self-activating gene regulatory circuit. The gene product X\ is its 
own transcription factor which, upon dimerisation, binds the promoter region of the gene thus triggering 
gene transcription. The transition rates corresponding to this gene regulatory circuit are given in Table |TV] 
For simplicity, we use an effective model in which the formation of the dimer and binding to the promoter 
region is taken into account in a single reaction, and the resulting number of promoter sites bound by two 
transcription factors is denoted Xi- 


Variable Description 

X] Number of transcription factor molecules 
X 2 Number of bound promoter sites in the gene promoter region 
X :i Number of unoccupied (unbound) binding sites in the gene promoter region 
Transition rate r Event 

W\(x) = R + k\X -2 r\ = (1,0,0) Synthesis of the transcription factor 

W 2 (x) = k‘>X\ ?’ 2 = (—1,0,0) Degradation of the transcription factor 

W 3 (x) = k 3 X\ (X\ — 1) V;> r .3 = (—2, +1, —1) Dimer binding to the gene promoter region 

W 4 (x) = A' 4 V '2 ?’4 = (+2, —1, +1) Unbinding from the gene promoter region 

TABLE IV. Random variables and transition rates associated to the stochastic dynamics of an auto-activation 
gene regulatory circuit 30 - 74 . X 2 corresponds to the number of transcription-factor dimer/promoter binding 
site trimers. See Fig. [6]for an schematic representation. 

cally represented in Fig. [6j In this circuit the gene product binds to form dimers which then act 
as its own transcription factor by binding to the promoter region of the gene. The rate-limiting 
factor is therefore the number of available binding sites within the promoter of the gene. For 
simplicity, our stochastic model associated to the rates shown in Table [IV] does not explicitly ac¬ 
count for dimer formation. We will assume that this process is very fast so it can be subsumed 
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under the formation of transcription-factor dimer/promoter binding site trimers (reaction 3, Table 
m. Furthermore, it is important to note that our stochastic dynamics exhibits a conservation law: 
X 2 + X 3 = e 0 at all time. This conservation law expresses the fact that the total number of binding 
sites, e 0 , is constant. 

In order to proceed with our analysis of the stochastic model of self-activated gene regulation 
(see Table [TV] and Fig. [6]), we apply the general methodology associated to our SCQSS approxi¬ 
mation. Following the general procedure explained in the previous sections, we start by deriving 
the stochastic Hamiltonian associated to the process defined by the transition rates shown in Table 
[TV] (see Section HD): 


H(p,Q ) = ( pi-l)(R + k 1 Q 2 p 2 ) + k 2 (l-pi)Q 1 + k 3 {p 2 -plp 3 )Q\Q 3 -\-k±(p 2 l p 3 -p 2 )Q 2 , (46) 

which, according to our theory (see Section HO), gives raise to the re-scaled Hamiltonian, H K (p, q ), 
defined by, 


H K (p, q ) = (pi - 1 )(R + n 1 q 2 p 2 ) + k 2 (1 - Pi)qi + {pi ~ P^P^q^-l'i + ^a{p\pz ~ ^ 2 )^ 2 , (47) 

where H(p, Q ) = k ?t ES 2 IIJj). q) and the re-scaled variables, q n and re-scaled rate constants, n v 
are defined in Table lYl 


Rescaled variables 

Dimensionless parameters 

r = k^ESt 

e = E/S, R = R/(k 3 ES 2 ) 

qi = Qi/S 

hi = ki/(k 3 S 2 ) 

Q2 = Q2/E 

k 2 = k 2 /(k 3 ES) 

© = Q3/E 

rc 4 = fc 4 / (k 3 S 2 ) 


TABLE V. Dimensionless variables used in Eqs. (l29l) . so is a characteristic scale associated to the aver¬ 
age number of molecules of transcription factor, X\ , and E is the average number of binding sites in the 
promoter of the self-activating gene. We further assume that S 2 > E. 

The re-scaled Hamilton equations are thus given by: 
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(48) 


e 


e 


e 


e 


dqi 

dr 

dq2 

dr 

dq 3 

dr 

dpi 

dr 

dp 2 

dr 

dpz 

dr 


R + hc 1 q 2 p 2 - «2<7 i - 2 qlq 3 p x p 3 + 2K 4 p 1 p 3 g 2 
(pi - i)ki ? 2 + q\q 3 - K A q 2 
- q{q 3 p\ + 

k 2 (Pi - 1) - 2<m 3 (P2 - P?P3) 

Kl(l - Pl)P2 - K4&3 - P 2 ) 

<?1 (PlP3 — P 2 ) 


(49) 

(50) 

(51) 

(52) 

(53) 


From these equations, we observe that for g 2 (r) + g 3 (r) = p, where p — e 0 /E, to hold we must 
have that p\{r) = 1 for all r. Imposing this condition on Eq. (1511) implies that p 2 (r) = p 3 (r), 
which, in turn, together with Eqs. (l52l) and (l53l) . imply that p 2 = p 3 =const. Finally, applying the 
QSS approximation to remaining equations, Eqs. (I48l)-(l50l). we obtain: 


dqi 

dr 


R + Ki pp 2 


dl 


+ q{ 


2 - «29l, 


qi=p-q 3 = p 




k 4 + q{ 


2 • 


(54) 

(55) 


As for the bistable enzyme-catalysed system, the parameter values are determined by matching the 
mean-field limit of our stochastic model, which is obtained by setting p % — 1 for all and p = 1 
(i.e. the number of binding sites exactly equal to its average), to the mean-field system proposed 
by Frigola et al.—. The mapping of our parameters to those of reference^ and their associated 
values are given in Table ED 

Finally, according to the theory developed in Section HU p 2 is determined in terms of the total 
number of binding sites within the gene promoter, eo'. 


dS 0 

-V 2 - 7 — = e 0l 
dp 2 


(56) 


where S 0 (p) = ln(G 0 (p)) and G 0 (p) is the generating function associated to the probability dis¬ 
tribution of eo, P(e 0 ). This probability distribution can be interpreted as corresponding to the 
distribution over a cell population of the number of binding sites in the promoter of gene x 4 . For 
example, if P(e 0 ) is a Poisson distribution the Eq. (l56l) reads^ 
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Rescaled parameter 

Parameter 

Units 

Reference 


K d = 10 

nM 

30 

K 2 = 1 

kdeg = 2 

min -1 

30 

K4 = 1 

r = 0.4 

nM • min -1 

30 

E> _ r 

kdegVTQ 

S = 1.0 



k^ES = kd e g 

E = 0.1 




TABLE VI. Parameter values used in simulations of the stochastic self-activation gene regulatory circuit. 



(57) 


where E = (e 0 ), i.e. the average of e 0 over a population of cells. Therefore, according to this 
analysis, we have that p = p 2 , provided that P(eo) is a Poisson distribution with parameter E. 


A. Bifurcation analysis 

Fig. m shows results regarding how the bifurcation diagram varies as we change pp 2 = pi, 
which, we recall, is determined by the (probability distribution of the) total number of binding 
sites within the gene promoter. Inspection of Eq. (l54l) shows that p 2 has the effect of renormalising 
the self-activation rate If p\ < 1 then the rate of gene self-activation is effectively reduce and, 
consequently the stability region of the inactive steady-state, qi ~ 0, is extended. That is, we need 
to go to larger values of k 4 to enter the region where the active steady-state, tp > 1 , becomes 
stable (see Fig. [7]). On the contrary, pi > 1 has the effect of extending the stability region of the 
active steady-state, q 1 > 1 . 

In order to verify the predictions of our bifurcation analysis (Fig. [7]), we consider Eqs. (l56l) 
and (l57l) . which relate the momentum variable p 2 to the number of binding sites within the gene 
promoter. If we assume that the latter is distributed according to a Poisson distribution, then Eq. 
(1571) holds and p 2 = p = e§/E. Under these conditions, our bifurcation analysis predicts that the 
probability distribution of X\, i.e. the random variable associated to the generalised coordinate 
qi, should change, as eo decreases, from being uni-modal with a single maximum about the ON 
value of Xi (or, when scaled with s 0 , qi) to exhibiting bi-modality, as the system approaches the 
saddle-node bifurcation which annihilates the ON state as it collides with the saddle point, with 
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FIG. 7. Bifurcation analysis for the SCQSS approximation of the stochastic auto-activation gene regulatory 
circuit Eqs. (I54l)- (l55l) . This figure shows the bifurcation diagram for different values of the parameters of 
P 2 - In these panels solid lines correspond to p 2 = 1, dashed lines to p\ = 0.9, dotted lines to p\ = 0.8, and 
dash-dotted lines to p\ = 0.7 (recall that p 2 = p). Parameter values as given in Table I VII 


two peaks about the ON and OFF states. If e 0 is further reduced the system will be driven passed 
this saddle-node bifurcation, the probability distribution becomes uni-modal but, unlike its large e 0 
counterpart, its peak is about the OFF -steady-state. We have verified this prediction by running 
simulations using the SSA. The results, which agree with our prediction, are shown in Fig. [8} 

Quantitative comparison between our asymptotic analysis and the simulation results follows 
the same procedure as in Section HTU i.e. we look at how the variance aforementioned behaviour 
regarding unbounded increase of fluctuations close to a bifurcation^ is used to locate the critical 
the variance a 2 = ((xi — (xi)) 2 ), with x\ = X\{r = T)/S changes as the control parameter 
varies: the maximum of a 2 as a function of the control parameter corresponds to the critical value. 
According to Fig. Ob), the critical value of p, p B , is approximately given by p B ~ 0.7. Our 
asymptotic analysis (see Fig. Ob)) predicts that p B = 0.78. 
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x/p 


FIG. 8. Simulation results for the stochastic gene regulatory circuit of self-activation (Table IfVl). We have 
plotted the probability P(x±,T ) = Prob(xi(r = T)) where x\ = X\/S and T = 100 as the number of 
binding sites in the gene promoter, given by X${t = 0) = pE. Average is performed over 1000 realisations. 
Parameter values are inferred from those given by Frigola et al.— as shown in Tables IVl and 1 VII We see the 
emergence of bistability at p = 0.7, whereas for smaller(larger) values of p, the system will be in the stable 
steady state corresponding to low(high) number of transcription factor molecules. 


V. CONCLUSIONS & DISCUSSION 

By means of the semi-classical quasi-steady state approximation, Section UH we have analysed 
stochastic effects affecting the onset of bi-stability in cell regulatory systems. Our theory shows 
that there exists a conserved momentum coordinate associated to each conserved chemical species. 
In the case of the enzyme-catalysed bistable system, Section HHl there are three such conserved 
momenta, associated to each of the conserved chemical species, i.e. Cdhl and its activating and 
inhibiting enzymes. For the self-activation gene regulatory network, we have one conserved mo¬ 
mentum, corresponding to conservation of the number of binding sites of the gene’s promoter 
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FIG. 9. Plot (a): Bifurcation analysis for the SCQSS approximation of the stochastic auto-activation gene 
regulatory circuit Eqs. (|54|)-(|55T). with K] = 3.0 Parameter values as given in Table I VII Plot (b): Simulation 
results for the variance a 2 = ((xi — (xi)) 2 ), with xi = X\(t = T)/S, associated to the simulation results 
shown in Fig. [U This plot shows how a 2 changes as the control parameter, p. The maximum of a 2 as a 
function of the control parameter helps us to quantitatively determine the corresponding critical value 73 . 


region. 

According to the SCQSS A analysis of22, the maximum rate achieved by an enzymatic reaction, 
V max , predicted by the mean-field theory^ is renormalised by a factor which equals the value of 
the (constant) momentum coordinate p, associated to the conserved enzyme: V™? = p ej PiVmax 
where Vmax is the maximum rate predicted by the SCQSSA. Similarly, we have shown that 
the mean-field maximum activation rate associated to the auto-activation gene regulatory model, 
A max , is renormalised in the presence of noise by a factor equal to the conserved momentum 
coordinate corresponding to the number of binding sites in the gene promoter, p 2 , i.e. Amax = 
pp 2 A max , with Amax being the SCQSSA maximum activation rate. As a consequence of this 
parameter renormalisation, we have shown that variation in the value of the conserved momenta 
can trigger bifurcations leading to the onset of bistable behaviour beyond the predictions of the 
mean-field limit, i.e. for values of parameters where the mean-field limit predicts the system to be 
mono-stable, the SCQSSA predicts bi-stability, and vice versa (see Figs. [2] and [7]) . 

Furthermore, we have established that the value of the constant momenta is actually determined 
by the probability distribution of the associated conserved chemical species, and, ultimately, by 
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the number of molecules of these species (see Eqs. (l44l) and (l56l)-(l57l)h Therefore, our theory 
establishes that the numbers of molecules of the conserved species are order parameters whose 
variation should trigger (or cancel) bistable behaviour in the associated systems. This prediction is 
fully confirmed by direct numerical simulation using the stochastic simulation algorithm (see Figs. 
[3l|4l and [8]). Quantitative comparison between the predictions of our asymptotic analysis and the 
simulation results (see Fig. [5] and [9]) shows that our theoretical approach slightly underestimates 
the critical value for the bistable enzyme-regulated system. The theoretical prediction for the 
self-activating gene regulatory network appears to slightly overestimate the critical value. 

Our results allow us to propose a means of controlling cell function. For example, regarding 
the enzyme-catalysed bistable model analysed in Section [TTT1, varying the number of molecules of 
the three conserved chemical species (Cdhl and the associated activating and inhibiting enzymes) 
enables us to lock the system into either of the Gi or the S-G 2 -M stable fixed points or to drive the 
system into its bistable regime where random fluctuations will trigger switching between these two 
states. This could be accomplished by ectopically increasing the synthesis of the corresponding 
molecule or by targeting the enzymes with enzyme-targeted drug s^— . Similarly, the dynamics 
of the self-activating gene regulatory system could be driven into or out of its bistable regime by 
supplying an inhibitor that irreversibly binds to the promoter region, thus decreasing the effective 
number of binding sites. 

This result allows us to explore strategies, for example, in the field of combination therapies 
in cancer treatment. Cellular quiescence is a major factor in resistance to unspecific therapies, 
such as chemo- and radio-therapy, which target proliferating cells. Bi-stability is central to control 
cell-cycle progression and to regulate the exit from quiescence, with enzyme catalysis (usually 
accounted for by (mean-field) Michaelis-Menten, quasi-steady state dynamics) being ubiquitously 
involve d-—™— . Our findings will allow us to formulate combination strategies in which chemo- 
or radio-therapy are combined with a strategy aimed at driving cancer cells into proliferation or 
quiescence depending on the phase of the treatment cycle. Evaluation of the viability and effi¬ 
ciency of such combination requires the formulation of multi-scale models^ 1 whose analysis is 
beyond this scope of this paper, and it is therefore postponed for future work. 

Our approach differs from previous work, such as Dykman et al.— in a significant aspect, 
namely, whilst their aim is to estimate the rate of noise-induced transition between metastable 
states in systems exhibiting multi-stability, the purpose of our analysis is to ascertain whether 
noise can alter the multi-stability status of the system. Dykman et al.— do not address such issue. 
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Eqs. (I36l)-(l39l) and (l54l)-(l55l) are derived from a semi-classical approximation of the Master 
Equation (or its equivalent description in terms or the generating function PDE). This approxi¬ 
mation yields a set Hamilton equations (Eqs. (8)-(9)) whose solutions are the optimal fluctuation 
paths and, as such, they describe fluctuation-induced phenomena which cannot be accounted for 
by the mean-field approximation. One of the best known examples of this is exit problems from 
meta-stable states in noisy systems (e.g. extinctions), where the semi-classical approximation pro¬ 
vides the optimal escape path from which information such as mean-first passage time or waiting 
time for extinction can be obtained (see, for example, references'^). Furthermore, Eqs. (1361) - 
(l39l) and (l54l)-(l55l) are derived from the general Hamilton equations, Eqs. (8)-(9), by means of an 
approximation based on separation of time scales, not on any mean-field assumption. 

A closely related subject to that analysed in this paper is that of noise-induced bifurcations^. 
Such phenomenon has been studied in biological systems where the mean-field limit does not 
predict bistability, such as the so-called enzymatic futile cycles^- where noise associated to the 
number of enzymes induce bistability. In the absence of this source of noise, the system does not 
exhibit bistable behaviour. We have not dealt with such noise-induced phenomena in the present 
paper, in the sense that all the systems analysed in this paper are such that their mean-field limit 
exhibits bistability. We leave the interesting issue of whether our SCQSSA framework can be used 
to analyse noise-induced bifurcation phenomena for future research. 
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